Thermally driven classical Heisenberg model in one dimension 
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We study thermal transport in a classical one-dimensional Heisenberg model employing a discrete 
time odd even precessional update scheme. This dynamics equilibrates a spin chain for any arbitrary 
temperature and finite value of the integration time step At. We rigorously show that in presence 
of driving the system attains local thermal equilibrium which is a strict requirement of Fourier law. 
In the thermodynamic limit heat current for such a system obeys Fourier law for all temperatures, 
as has been recently shown [A. V. Savin, G. P. Tsironis, and X. Zotos, Phys. Rev. B 72, 140402(R) 
(2005)]. Finite systems, however, show an apparent ballistic transport which crosses over to a 
diffusive one as the system size is increased. We provide exact results for current and energy profiles 
in zero- and infinite-temperature limits. 

PACS numbers: 44.10.+i,75.10.Jm, 66.70.Hk 



I. INTRODUCTION 

The flow of heat from a hot source to a cold sink is 
conventionally described in the hydrodynamic limit by 
Fourier law J = -kVT, where J is the steady state 
thermal current set up in response to the temperature 
gradient VT and k is the (finite) thermal conductiv- 
ity. There have been several attempts [HQ at a mi- 
croscopic 'derivation' of this phenomenological equation 
which, however, could not be achieved yet. Needless to 
say, in spite of the huge amount of studies over decades, 
our understanding of this basic transport phenomenon is 
still not quite satisfactory. It has been found that in a 
variety of one-dimensional models [Tj-Q the thermal cur- 
rent scales with system size as where a < 1. This 
corresponds to a diverging heat conductivity in the ther- 
modynamic limit and thus is a violation of Fourier law. 
It is quite a surprise that Fourier law, which has been re- 
markably consistent with experimental results in general, 
is found to be invalid in many models in low dimension. 
For three-dimensional systems Fourier law is believed to 
be generically true but a rigorous proof is still lacking Q . 

In a recent work it has been shown numerically (us- 
ing the Green Kubo approach Q) that a classical one- 
dimensional Heisenberg spin model obeys Fourier law 
at all temperatures yjj. Also, it has been known for 
quite some time now that at infinite temperature such 
spin systems follow the energy diffusion phenomenology 
in the hydrodynamic limit On the other hand, the 
ID spin-i quantum Heisenberg model (QHM) being in- 
tcgrable, violates Fourier law and thermal transport is 
ballistic 0, 13 ■ For recent reviews on the theoretical and 
experimental developments in quantum spin models see 
@ and references therein. One of the basic assump- 
tions for the validity of Fourier law lies in the estab- 
lishment of local thermal equilibrium (LTE) in the sys- 
tem 0, H, [Io| , which allows one to define thermodynamic 
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quantities in the steady state such as pressure, temper- 
ature etc., locally. Most studies, however, focus on the 
issue of whether Fourier law is obeyed or not, without ex- 
plicitly verifying the existence of LTE. In fact, it is known 
that a system may settle down to a noncquilibrium steady 
state (NESS) which does not satisfy the essential require- 
ment of having LTE e.g., XY model, Lorentz gas model 
[lH ] . It is believed that the absence of LTE in these ex- 
amples is due to the existence of infinitely many local 
conserved quantities in the dynamics [ll|. Many of the 
theoretical approaches @, Q rely on Linear response the- 
ory (Green-Kubo formula) where conductivity is mea- 
sured by computing two point current-current time cor- 
relation which assumes quasi-equilibrium. The Kubo for- 
mula [H is strictly valid close to equilibrium and in the 
limit L — > oo, and so considerable care should be taken 
in making conclusions from experimental or simulation 
data which deal with finite system size and drive Q. 

In this paper, we study thermal transport properties of 
a classical one-dimensional Heisenberg spin model. We 
use a discrete time odd even (DTOE) dynamics which, 
unlike standard numerical integration schemes, evolves 
the system to the correct steady state without violating 
the required conservations. We explicitly show that the 
DTOE dynamics equilibrates a closed system and the fi- 
nal state is the same for all nonzero values of At. With 
two equal temperature baths attached to its two ends, 
the system eventually equilibrates under the DTOE dy- 
namics and attains the temperature of the baths. When 
temperature of the heat baths is different, thermal equi- 
librium is established in the system locally. With finite 
drive, we study in details the transport properties of the 
system e.g., thermal current J, energy profiles, conduc- 
tivity k without invoking linear response theory. We 
find that, in the thermodynamic limit, the system obeys 
Fourier law at all temperatures which is consistent with 
the recent study For a finite system however, there 
is a characteristic temperature below which the system 
crosses over to a regime where transport becomes bal- 
listic. We present exact results for thermal current and 
energy profile in the limit T — > and T — > oo. 
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The paper is organized as follows. In Sec. |TT] and Sec. 
Ill Al we describe the model and the spin dynamics in de- 
tail. We then look into the equilibration of a closed spin 
chain under this dynamics in Sec. IIII1 In Sec. IIVI we 
study, analytically and numerically, the transport prop- 
erties of the model in presence of thermal baths. We 
present a discussion and summarize our main results in 
Sec. H 



II. MODEL 

Consider classical Heisenbcrg spins {Si} (three- 
dimensional unit vectors) on a one-dimensional regular 
lattice of length L (1 < i < L) with periodic boundary 
conditions. The microscopic Hamiltonian is given by 



L 



'i+l 



L 

-K cos i 



(1) 



where the spin-spin interaction is ferromagnetic for cou- 
pling K > and anti-ferromagnetic for K < 0. For all 
the numerical results shown in the paper K has been set 
to unity. Here Qi is the relative angle between Si and 
Si+\. The microscopic equation of motion can be taken 
as 



S> x Bi 



(2) 



where Bi = K(Si^\ + Sj+i) is the local molecular field 
experienced by the spin at site i. Clearly, Eq. @ con- 
serves (i) the magnitude of the individual spin vectors Sf 
and (ii) the energy density 



section, we show that a straightforward numerical inte- 
gration of Eq. ([2]) fails to conserve either Sf or E or 
both. We also discuss in detail the advantages of using 
the discrete time odd even dynamics (DTOE). 



A. Why DTOE dynamics? 

To integrate the equation of motion numerically one 
would naively start off with a finite difference equation 
of the form 



S. 



i,t+l 



Si. 



At 



SxB 



(4) 



and update all the spins at time t to obtain their values at 
i+l. Such an Eulerian scheme cannot be used for this 
system because of the fact that it does not satisfy the 
required Sf and E conservations. It is easy to calculate 
the energy E(t) and Sf(t) in the Euler scheme using Eq. 
which comes out to be 



t-1 L 



E(t) = E(0)-K(At) 2 Y / J2[ SxS 

T=0 1=1 
t-l 

S?(t) = S?(0) + (AtfY,[SxB 



SxB 



i+l,T 



T=0 



2 

y.r 



(5) 



Thus, any finite At, however small, breaks both the con- 
servations and consequently the scheme fails for all prac- 
tical purposes [see Fig. [Ha)]. A way to keep the mag- 
nitude of the spin vectors conserved is to use a spin pre- 
cession dynamics 



E 



1 L 



where e, 



-KSi ■ Si+i. 



(3) 



Note that this dynamics is the classical equivalent of the 
quantum dynamics for a spin-i QHM. Just as the com- 
mutation relations of quantum spin operators, the classi- 
cal spins components obey the standard Poisson bracket 
relations for angular momentum. 

However, there is a fundamental difference between the 
quantum spin-^ model and the classical model. The spin- 
\ QHM is integrable, whereas all higher spin [S > 1) 
QHMs (and therefore the classical model which corre- 
sponds to S — > oo) are non- integrable. Consequently, 
there are infinitely many conserved quantities in a spin- i 
QHM (which includes the energy current) and the ther- 
mal transport is ballistic 0]. On the other hand, only the 
total spin and the total energy are conserved in the corre- 
sponding classical model, and thus one expects transport 
properties to be normal. 

Since we wish to study thermal transport, typically far 
away from equilibrium (for which no general theoretical 
formulation is known), we need to integrate Eq. @ nu- 
merically keeping the conservations intact. In the next 



Si. 



S cos 4> + (S x B) sin <j> + (S- B)B{1 - cos < 



(6) 



where Bi = and <fo = \Bi\At [1J], instead of 

Eq. (|2|). A spin Si, when updated using Eq. under- 
goes a precessional motion about the instantaneous local 
molecular field Bi which keeps its magnitude unaltered 



Sf t . However, this precessional dynamics 



does not preserve energy conservation. Expanding Eq. 
© in powers of At and retaining terms up to O(At) one 
obtains back the first equation of Eq. ([5]) and thus energy 
conservation still remains violated for any Ai > 0. This 
has also been shown numerically in Fig. QJb). Other 
numerical schemes such as Rungc-Kutta, predictor cor- 
rector method etc., will also fail to preserve the energy 
conservation for the same reason. 

We now describe an odd-even update rule which along 
with the precessional dynamics has been herein referred 
to as the DTOE dynamics. Starting from a spin con- 
figuration {Si}, we numerically implement the dynamics 
described in Eq. (1) by alternate parallel updates of the 
spins on odd and even sublattices. Henceforth, we refer 
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to these two groups of spins as odd and even spins. At 
each Monte Carlo step (MCS), first only even spins are 
updated using the spin precession dynamics Eq. © and 
the odd spins are kept unaltered. Next, the spins on the 
odd sublattice are updated similarly. These two steps up- 
date all the spins {Si} in the system and constitute one 
MCS. It is straightforward to check that update of any 
spin Si affects only the energy of the neighboring bonds 
£i_i and ej but their sum (ej_i + e,-) remains constant. 
Thus DTOE dynamics is strictly energy conserving. 



Si 2 (t) 




E(t) 






(a) 




2500 5000 7500 10000 
t 



2500 5000 7500 10000 
t 



1.03 


S, 2 (l) 


-0.42 


1.02 


E(t) y / 








-0.43 


tn" 1.01 


(C) 




1 


-0.44 


0.99 




-0.45 



1.03 


Si 2 (t) 




-0.42 


1.02 


E(t) 




-0.43 


1.01 




(d) 




1 




-0.44 








0.99 




-0.45 



2500 5000 7500 10000 
t 



2500 5000 7500 10000 
t 



FIG. 1. (Color online) Evolution of spin magnitude Sf and 
the energy density E. (a) Parallel update using Eq. Q (b) 
Parallel update using Eq. ((6} (c) Odd even update using Eq. 
© (d) DTOE. For all the figures At = 0.001 and L = 10000. 
Thus, only DTOE dynamics conserves both E and Sf . 

Clearly, the spin precession dynamics conserves the 
magnitude of the spin vectors while, energy conservation 
is maintained by the odd-even update rule (see Fig. [TJ). 
A recent paper [14[ has also employed this odd-even pro- 
cessional dynamics with large At to study transport in a 
classical Heisenberg model (ID periodic spin system) in 
presence of quenched disorder numerically. Although this 
dynamics does not directly follow from the equation of 
motion [Eq. @], it can be used to study the system nu- 
merically provided that the system equilibrates for any 
arbitrary At. In the next section, we study the equi- 
libration of a closed system when evolved using DTOE 
dynamics. 



correlation functions are obtained from a closed system 
with a fixed energy (i.e. in a micro-canonical ensemble). 

The partition function of the system [l3| with the 
Hamiltonian given by Eq. @ is 



Z 




(7) 



i=l 



j3 = l/ksT and ks has been set equal to unity hence- 
forth. The two-spin correlation functions are therefore 
given by 



C lr = (Pl(Si- Si+r)) 



HflK) 



where Pi are Legendre polynomials and 



(8) 



Pl(x)e kx dx. 



(9) 



These correlation functions can be written explicitly in 
terms of the Langevin function C(x), for example, 



Clr = WK)\ 



C 2r = [1 - 3C(I3K)/((3K )} r . 



(10) 

It is evident that average energy of the system is —KCn 
and thus, a closed system with a fixed energy density E 
has an effective /3 = Cr x {-ElK)jK. 

First, let us compute the equal-time spin-spin corre- 
lations Ci r (t) for a closed system and check that these 
evolve to the stationary value given by Eq. (|10[) . The 
time series for C\ r (t) and C2 r (t) with different r values 
are shown in the Fig. [2j We find that the correlation 
functions for different At saturate to the same value at 
late times. In Fig. [3] the time averaged equilibrium corre- 
lation functions C\ r obtained from systems set at different 
energies, are also found to be in remarkable agreement 
with Eq. (|TU|) . Thus the DTOE dynamics equilibrates 
the system irrespective of the value of At; At only al- 
ters the equilibration time of the system. A larger At is 
preferable as equilibration in this case is attained faster. 



IV. OPEN SYSTEM 



A. Modeling heat bath 



III. CLOSED SYSTEM 

We first investigate whether a closed system (i.e., with 
periodic boundary conditions) under DTOE dynamics 
evolves to the correct steady state for different values 
of At. To do this, first we compute the correlation func- 
tions of the system in canonical ensemble, subjected to 
temperature T and then show numerically that the same 



In order to study energy transport in the system, wc 
now look into an open system with heat baths attached 
to its two ends. The left and right baths are set at tem- 
peratures 1 1 Pi and 1 / ' f3 r respectively. Each bath is known 
as a stochastic thermal bath [l5j which means that it is in 
equilibrium at its respective temperature and has a Boltz- 
mann energy distribution. The baths are implemented by 
introducing two additional sites i = and i = L + 1 in 
the system with spins Sq and Sx+i respectively. These 
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FIG. 2. (Color online) Semi-log plots for the evolution 
of equal-time correlations Cn, C12, C13, C21, C22, C23 
as obtained from simulation for a closed system of L = 
1000. DTOE dynamics has been used with time-step At = 
1.5, 1.0, 0.75, 0.5, 0.25 and E = -0.5. 

pairs of spins (5*0, Si) and (Sx, Sl+i) behave as stochas- 
tic heat baths at two ends of the system. The baths are 
in equilibrium at their respective temperatures and the 
bond energies eo and £l have a Boltzmann distribution 

P(e ) ~ e" Aeo and P(e L ) ~ e"^ . (11) 

The interaction strength of the bath spins with the sys- 
tem is taken to be A, and therefore both eo and e^ are 
bounded in the range (—K,K). Thus the mean energies 
of the left and the right bath are given by 

E t = (e Q ) = -K£(ftK) and E r = (e J = -KC{p r K). 

Following the odd-even rule, the spin So is updated 
along with the even spins, whereas Sl+i is updated with 
odd (even) spins depending on whether L is even (odd). 
To update So, first the energy of the bond eo between 
the spins (So, Si) is set to a value drawn randomly from 
P(eo) given in Eq. (jlip . The spin So is then constructed 
such that eo = —A' So -Si. During this update Si is not 
modified as it belongs to the odd sublattice. At the right 
end, the spin Sl+i is updated similarly. We must men- 
tion that energy conservation is violated during update 
of bath spins. The interaction of the bath spins with the 
neighboring spins allow boundary fluctuations to prop- 
agate into the bulk, thus inducing a thermal current in 
the system. 
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FIG. 3. (Color online) Semi-log plot for different equilibrium 
correlations Cw for I = 1, 2 vs. C11 = —E obtained from 
simulation (points) and compared with theory (lines) for a 
periodic lattice of L — 1000. DTOE dynamics has been used 
with time-step At — 1.0. 

Although the closed system equilibrates under DTOE 
dynamics, it is not guaranteed that an open system will 
also equilibrate when baths are attached to its two ends. 
Before studying the system with a finite drive, we study 
the equilibration of an open system with baths main- 
tained at the same temperature, thus still keeping the 
system in equilibrium at temperature T, i.e. /3; = /3 = j3 r . 
With baths maintained at equal temperature T, the spin 
chain is expected to eventually reach a thermodynamic 
equilibrium corresponding to the bath temperature T. 
We have calculated numerically the average energy of 
the system (E) for different values of f3 — 10, 1, and 0.1, 
which is shown in Fig. 0Ja). Evidently, at late times (E) 
approaches the stationary value —KjC(K(3). Figure |U[b) 
shows that the system attains a unique stationary state 
consistent with the bath temperature and this final state 
is independent of the value of At used. Similar to the 
case of the closed system, At decides only the equilibra- 
tion time of the system. We also measure the equilibrium 
correlations C; r which are shown in Fig. @Ic) where the 
numerical values are found to be in agreement with Eq. 
(fTU|) . This assures us that for any nonzero At the DTOE 
dynamics is no different from the equation of motion (f2|) , 
and allows the system to attain the correct equilibrium 
state. 



B. LTE in a driven system 

A finite thermal drive is imposed on the spin chain 
by setting the two heat baths at unequal temperature, 
i.e. pi ^ j3 r . The bath and bulk spins are updated as 
mentioned in the previous section. Now, since the bath 
temperatures are unequal, the system is driven out of 
equilibrium. However, it may still be possible to define 
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FIG. 5. (Color online) Semi-log plot for different correlations 
Cir(x), where < x = i/L < 1, for I = 1,2 vs. Cn(x) ob- 
tained from simulation (points) and compared with Eq. (|10[) 
(lines) for a open system of L — 1000. The average bath en- 
ergies are Ei — —0.7 and E r — —0.1. DTOE dynamics has 
been used with time-step At = 5.0. 



FIG. 4. (Color online) Semi-log plot for the evolution of 
average energy {E) for a system with j3i = /3 r = /3 for 
(a) P = 0.1, 1.0, 10.0 and a fixed At = 1.0 and (b) At = 
1.0,0.5,0.1,0.05 and a fixed /3 = 1.0. (c) Different correla- 
tions Cw obtained from simulation (points) are shown as a 
function of Cn along with the functions Eq. (|10|l (lines); 
L = 1000 and At = 1.0. 

a temperature like thermodynamic variable locally in a 
region if the average energy of the sites (e,) belonging 
to that region is not too different from each other. The 
system is said to have local thermal equilibrium (LTE) 
if all the correlation functions measured in this local re- 
gion are identical to those of a thermodynamically large 
equilibrium system with an average energy (ej). 

Numerically, one measures the correlation functions 
Ci r {x) locally over n <C L consecutive sites about x — 
i/L such that the average energy of these n sites is almost 
equal to each other. For a system of size L = 1000 we 
measure Ci r (x) up to three nearest neighbors for / = 1, 2 
and by averaging them over n = 20 sites. This is shown 
in Fig. [3J where we have shown Ci r (x) for different x 
(in the range (0, 1)) as a parametric function of C\\{x). 
For comparison, the equilibrium curves (from Eq. (|10[1 ) 
are also shown in the figure as solid lines. An excellent 
match with the equilibrium functions assures that the 
driven spin system attains thermal equilibrium locally. 
Thus, following the equilibrium definition, we may de- 
fine uniquely the local inverse temperature 

p{x) = ±-£-\ Cl i{x)). (12) 
A 



C. Analytical Results 

In the previous section we have seen that the driven 
system attains local thermal equilibrium, and thus a local 
temperature can be uniquely defined through Eq. (|12|) . 
Therefore, the usual definition of Fourier law J oc VT can 
also be equivalent ly expressed as J oc VA. The thermal 
current and the energy profile are measured as described 
in the following. Since the DTOE dynamics alternately 
updates only half of the spins (but all the bond energies 
simultaneously), the energy of the bonds e° measured 
immediately after the update of odd spins is different 
from ef measured after the update of even spins. Clearly, 
this difference t\ — e° is a measure of the energy flowing 
through the i-th bond in each MCS. Thus the thermal 
current in the steady state is given by 

J = {el - el) (13) 

and the average energy of i-th bond is = \ {e\ + e°). In 
fact, this expression for the current, in the limit At — > 0, 
is consistent with that obtained from the continuity equa- 
tion ii{t) — Ji-i(t) — Ji(t) where = —KSi ■ Si+\ is 
the local energy density. Straightforward calculation us- 
ing the continuity equation gives the instantaneous cur- 
rent across z-th bond 

Ji(t) = K Si ■ (S i+1 x S i+2 ). (14) 

Again, when the i-th site (say, even) gets updated, the 
energy of the i-th bond is ef = —A Si(t + 1) • Si + i(t) 
and after the subsequent update of odd sites it becomes 
e° = -A Si(t + 1) • S i+1 (t + 1). Thus using Eq. ©, the 
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instantaneous current (time is measured in units of At) 
across the z-th bond in the limit At — »■ 0, reduces to 

J t (t) = el - e° = K S t ■ (S i+ i x S i+2 ), (15) 

which is same as Eq. (|14l) . 

In the following, we show that the model with DTOE 
dynamics can be solved exactly in both T — s> and 
T — > oo limits to obtain analytical expressions for the 
energy current and energy profile. We show that en- 
ergy transport in any finite system is ballistic in the 
limit T — > 0, whereas diffusive transport is observed for 
T -> oo. 

Ballistic Limit (T — > 0): In this limit, the spins are 
nearly aligned (hence Si x Bi ~ 0) and, therefore, pre- 
cess by small angles. In DTOE dynamics, the angle of 
precession is </>,; = |_B, |At and so for small At it mimics 
the dynamics of the system at low temperature. This 
equivalence can be utilised to write the energy function 
Eq. in this limit as 

H^-Jf^l-^ 2 ), (16) 

which is similar to the energy function of a harmonic 
system. Since the DTOE dynamics is energy conserving, 
update of a spin Si assures that 0f_ 1 + Of remains invari- 
ant. Again, since the low temperature stationary state 
dynamics is governed by spin waves i.e., spin configura- 
tion whose orientation varies slowly with distance along 
the axis of the chain, the spins are locally parallel and 
Oi-i + 8i is invariant in the At — > limit. Thus the only 
allowed dynamics for the angle variable is 

Qi-i,t — > &iJ+At &i,t — » Oi-i,t+At- (17) 

Consequently, the energy of the bonds that connect to the 
i-th spin, namely e%-\ and ej, are mutually exchanged in 
this limit. Starting from t = 0, the bond energies £2i+i 
for the odd bonds 'move' to the right and the even bond 
energies e2i 'move' to the left ballistically (without any 
scattering) . 

In the steady state, the average bond energies after 
the update of odd spins are (e^) = Ei and (e% i+1 ) = E r 
whereas, the same after the update of even spins become 
(e|j) = E r and (e| i+1 ) = Ei. Thus current J = E\ - E r is 
independent of the system size L and thermal transport 
is ballistic. In the steady state, the temperature is same 
at all bulk sites, which is given by 

Tbulk = jr-i ( E,+E7 j ( 18 ) 

Note that this ballistic behavior is a consequence of the 
fact that the limit At — > is taken before the thermo- 
dynamic limit L — > oo. When At ~ 0, the spins pre- 
cess slowly since the precision angle <pi is proportional to 
if At. Thus the effective correlation length diverges as 



At — > and energy in any finite system would be trans- 
ferred to arbitrary distances without being scattered. 

Also, unlike the equilibrium case, the correlation func- 
tions of the driven system depend on At when L is fi- 
nite. In the following, we argue that when At — > the 
correlation length £ actually diverges. Since the energy 
profile in this limit is flat, a small change At —5- At' will 
not change the correlation functions substantially. In or- 
der to keep the correlation functions (which are functions 
/3K) unaltered, the inverse temperature /3 should scale as 
/?' = (3At/Af, so that K'fi = K/3. Again, since the cor- 
relation length £ = | \ n c{Kp)\ (calculated from Eq. (flO|) 
taking Ci r = e~ r ^) diverges linearly with /3 in the limit 
P — > oo, we have 

e-(At)- 1 . (19) 

This indicates that the steady state energy profile also 
depends on At, which will be discussed later in section 
HVlE (see Fig. [TUJa)). We must mention here that this At 
dependence is only a numerical artifact in finite systems. 
In fact, 1/At has to be compared with the two other 
length scales of the problem, namely, the size of the sys- 
tem L and 1/T (as the correlation length also diverges 
in T limit) and thus, the effective correlation length 
will appear to be 1/At only when both L and 1/T are 
much smaller. In other words, for the numerical integra- 
tion of Eq. , one must choose the integration time step 
At larger than both 1/L and T to avoid dependence of 
the steady state on At. Therefore, a thermodynamically 
large system in this problem corresponds to a system with 
L > 1/T > 1/At. In this limit, the correlation length £ 
remains smaller than L for all T > 0; the steady state be- 
haviour is independent of At and one recovers diffusive 
thermal transport (see Fig. [TT] and related discussions 
later) . 

Diffusive limit (T — > oo): In the other limit T —> oo, 
the spin orientation is random and thus the dynamics is 
equivalent to large At limit, where the precession angle 
</> is large and effectively the spin precesses by a random 
angle. Updating the i-th spin then results in random re- 
sharing of the bond energies e^-i and ej obeying the local 
energy conservation imposed by the DTOE dynamics. 
Effectively, 

64-1,4+1 = r(ei_l -I- €i) t £i,4+l = (1 - r)(ej_i + e.;) t 

(20) 

where, r is a uniform random number in the range (0, 1). 
This dynamics is similar to the diffusive dynamics dis- 
cussed by Kipnis et. al. [l(| except the fact that here we 
use the DTOE dynamics. In the steady state, the average 
energies at different sites satisfy the following equations. 
Update of odd sites ensures that for j = 0, 1, . . . , L/2, 

(4i) = «4> + /2 

<4,+i> = «4> + <4+i» /2 (2i) 
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Similarly update of odd sites gives 

= ((4i-i) + /2 
(4) = ((4i-i) + (<&» A (22) 

for j = 1, 2, ... , L/2, along with the boundary conditions 
(e§> = E l (e e L+1 ) = E r . (23) 



These set of linear equations (j2T ]) -([23 |) provide a unique 
solution 



(24) 



where r\ = 0, 1 for i = even, odd respectively. Clearly the 
energy profile e, = ((e° + e|})/2 is linear and the current 



J = 



E r — Ei 
L+2 



follows Fourier law. 



D. Thermal current 

For finite T, the model is not analytically solvable and 
we study transport properties numerically using DTOE 
dynamics. Two thermal baths are attached to the two 
ends of the system having average energy E\ = E and 
E r = E + AE respectively. The steady state current J, 
measured using Eq. ([T3]). is shown in Fig [6] Clearly, J 
decreases with increase of system size, and approaches 
the algebraic form J ~ 1/L in the thermodynamic limit. 
For small L, however, J varies slower than 1/L. Keeping 
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FIG. 6. (Color online) Log-log plot of the steady state current 
J vs. L agrees well with Fourier law J ~ L^ 1 for large L, but 
deviates for small L values. The average energy of the two 
baths are E t and E r = E\ +AE with AE = 0.1 and At = 1.0. 

this in mind, a suggestive phcnomenological equation for 
the energy current can be written as 



J 



AE 



(25) 



where n and £ are parameters which depend on T, At. As 
the temperature T — > 0, the correlation length of the spin 
chain £ — > oo, and consequently heat transport shows an 
apparent ballistic behavior. In the other limit, i.e. when 
T is large and £ — >■ 0, thermal transport in the system is 
diffusive. 

Following Eq. (|23|) , k and £ can be measured from the 
slope and intercept of the straight line L = K^- — £. 
In the inset of Fig. [71 we have shown L against AE/ J 
for different bath temperatures; all the curves are linear 
and the best fitted straight lines give respective k and £. 
Further, we observe that the parameters k and £ always 
maintain a fixed ratio with each other for any given At. 
This becomes evident from the collapse of J versus L/k 
curves for different bath temperatures (see Fig. [7J. This 
implies that k should have the same T dependence as £. 
Since near T = the correlation length £ ~ T -1 , we 
expect that k should also diverge inversely with T in the 
limit T — > . In fact, k is the conductivity of the system 
in the thermodynamic limit L £ and its divergence at 
T = indicates that the system is near a critical point. 

The behavior of k with temperature T is shown in Fig 
[5] Close to T = the system relaxes extremely slowly 
and numerical studies in this limit become computation- 
ally expensive. One needs to go to extremely small tem- 
peratures to see the T _1 divergence of k, which could not 
be reached with the available computational resources. 
The inset of Fig. [S] shows that k vanishes linearly in 
the limit At — >• 0, which can be understood as follows. 
The DTOE update of a spin 5, keeps the local energy 
(ej_i + e;)/2 conserved, i.e., de/dt = 0. Therefore, for 

finite At we must have Ae ~ (At) 2 so that lim = 0. 

Again from Eq. ([TB"]) we have 



At^o At 



J - (Ae) - (At)' 2 



(26) 



Since £ diverges as (At) 1 (from Eq. (fT9|) ) and the cur- 
rent in this limit J ~ j, we have n ~ At. Until now 
we have discussed thermal transport for a small AE and 
assigned the k, obtained from Eq. (1231) . to be the con- 
ductivity of the thermodynamic system at energy E (or 
temperature T). As such, in this limit the system is not 
too far from equilibrium, in a way that all parts of the 
system are maintained almost at the same temperature 
T. However if AE is appreciably larger, both local energy 
and its gradient varies significantly across the system. In 
such a case, one can appropriately define a local conduc- 
tivity as, 



Klo 



J 



de(x) 
dx 



(27) 



local 



To measure Ki oca i, we set the bath energies at Ei and 
E r -c Ei and calculate the energy profile e(x) and its 
de(x) 

gradient — ; — at different x = i/L along the system. 
dx 

The inset of Fig. [9] shows the energy profiles obtained 
for different bath energy E\ and E r = E\ — 0.5. In the 
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FIG. 7. (Color online) Collapse of the curves J vs L/k, for 
different values of average bath energies Ei and E r with a fixed 
AE = 0.1. The simulation data are shown by points and the 

AE 

solid line corresponds to the curve of the form , _ , . ti— ■ 

(L/k) + Z/k 

The k values are obtained from a straight line fit of the form 
L = n^S- — £ as shown in the inset (see text). 
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FIG. 8. (Color online) Variation of k with temperature T, 
as obtained from the straight line fit using Eq. (|25l) . shows a 
divergence in k as T — > 0. (Inset) k varies linearly with At. 



main figure, we have shown the local conductivity K,i oca i 
as a function of the local energy e(x); the overlapping 
regions, although obtained from energy profiles with dif- 
ferent boundary energies, match remarkably. Thus K[ oca i 
is a well defined function of the energy (or equivalently, 
temperature) and has the same value for a given energy, 
irrespective of the average energy of the two baths. Since 
the spin system attains local thermal equilibrium for all 
nonzero temperatures, Ki oca l at a given local energy e 
must be same as the conductivity k calculated using Eq. 
(|25p for a large system with average bath energies e and 
e + AE respectively. This is shown as open circles in Fig. 



Hfor AE = 0.1, and different e. 




-1 -0.8 -0.6 -0.4 -0.2 



FIG. 9. Variation of local conductivity Ki oca i with energy 
e(x). Ki oca i has been calculated locally (from the profiles 
shown in inset) using Eq. (|27[) . In the overlapping regions 
of the energy profiles, Ki oca i for different energies collapses 
onto a single curve. As energy becomes small (which corre- 
sponds to T — > 0), the local conductivity diverges. The open 
circles correspond to k calculated using Eq. (|25|) . AE = —0.1, 
and for different Ei = —0.2 to —0.9 in steps of —0.1. For both 
the figures At = 1.0. and L = 1000. 



E. Energy profiles 

We now turn to the energy profile of the driven system 
and investigate the dependence of the same for the 
following three cases: 

a. At dependence. The energy profile for a finite 
system depends on the parameter At as can be seen from 
Fig. [TW a). We have shown earlier (see section lTVl C) that 
for finite L, the two asymptotic limits At — > 0, oo corre- 
spond to ballistic and diffusive transport respectively and 
hence it is expected that for a smaller At the profile will 
be relatively flatter as compared to a larger value of At. 
Thus for any finite system if L <C £ the transport will 
be ballistic (fiat energy profile) and one has to simulate 
larger systems to observe a diffusive behavior (linear en- 
ergy profile). 

b. E dependence. A lower E implies a lower temper- 
ature T and hence for a given L and At, the correlation 
length monotonically increases as E is decreased. The 
system approaches a ballistic limit with energy profile as 
E is reduced for a given value of At and L. However in 
the thermodynamic limit and for T > 0, Fourier law is 
always satisfied. This is shown in Fig. Q3B» 

c. L dependence. The L dependence of the energy 
profile is also consistent with what we have already dis- 
cussed. For a given value of At and E, a smaller L shows 
a flatter profile as compared to a system with larger L, 



9 



as can be seen from Fig. HUTc). 
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x = HL 

FIG. 10. (Color online) Energy profiles for (a) different At 
with fixed bath energies Ei = —0.2, E r = —0.25 and L — 100. 
(b) different bath energies E u E r with AE = -0.05, At = 0.1 
and L — 100. The profiles have been shifted up along the 
energy axis by \Ei \ to accommodate all the profiles within the 
same energy window, (c) different L with fixed bath energies 
Ei = -0.6, AE = -0.1 and At = 0.1. 



V. DISCUSSION 

To summarize, we have studied thermal transport in 
a one-dimensional classical Hcisenberg spin model using 
discrete time parallel even-odd updates with spin pre- 
cession (DTOE) . While conventional integration schemes 
fail to preserve the required conservation of Sf and E, 
this dynamics preserves both. The DTOE dynamics con- 
verts the equation of motion fl^J) to a map (Eq. ([6])) 
with an additional parameter At (besides the interac- 
tion strength K and the system size L); the equation of 
motion is recovered from the map in the limit At — > 0. 
We explicitly show that this energy conserving dynam- 
ics equilibrates a closed system (having a fixed energy) 
and an open system attached to equal temperature heat 
baths, for any finite At. When the system is driven by 
maintaining a finite temperature difference between the 
two ends, we explicitly show that the system attains local 
thermal equilibrium. However, the steady state proper- 
ties such as the correlation length £ (Eq. (IT9|) ). thermal 
current (Eq. (|26])). and energy profile (Fig. [10] (a)) de- 
pend on At when the system size L is finite; such spurious 



At dependence disappears in the thermodynamic limit. 

Our numerical simulations of the system for different 
bath temperatures suggest that the thermal current J 
can be expressed in the form J = n-j^, where k and 
correlation length £ depend on the temperature T and 
At. In the thermodynamic limit L 3> £, the spin sys- 
tem exhibits Fourier law; the energy profile in this case 
is linear with the slope asymptotically approaching the 
value to* = Ap. However, for small system sizes (i.e. for 
L -C £) , thermal transport appears to be ballistic with a 
relatively flatter energy profile. The same scenario pre- 
vails when, instead, the temperature T is varied. That 
is, for a given L and At, the slope of the energy profile 
approaches to* (or 0) as T — > oo (or 0). Thus, finite 
systems show an apparent crossover from a diffusive to a 
ballistic behavior as T is lowered below a characteristic 
temperature scale T*. This is described in Fig. QT] along 
with additional numerical evidences. 

To demonstrate the crossover phenomena quantita- 
tively, we measure the local slope to of the energy profile 
at i = L/2 for a system of size L and fit it to a functional 
form 

m(T,At)=m* T /A ~ . (28) 

v > > T + T*(At) v ' 

Clearly, T* is the value of temperature for which the slope 
is half the desired slope for diffusive transport, m* = Aj^. 
In Fig. [TT]we show the crossover temperature T*(At) for 
two different system sizes, L = 100, 200. The crossover 
line T*(At) separates the diffusive regime (well above 
the curve) from the ballistic one (well below the curve). 
Again, the crossover line shifts downwards when the sys- 
tem size L is increased. This clearly indicates that in the 
thermodynamic limit, the crossover line is infinitesimally 
close to the axes and thus, for any T > 0, one observes a 
diffusive behaviour independent of the choice of At. The 
apparent ballistic behaviour in the small At or small T 
regime is only an artifact of finiteness of the system and is 
a consequence of the divergence of the correlation length, 
as £ ~ 1/T and £ ~ 1/ At. The dependence of £ on At 
can be understood from the fact that the spin precesses 
by a small angle, proportional to At. This effect is simi- 
lar to a low temperature precession-dynamics, where the 
correlation length £ is very large. Thus, for studying 
thermal transport at a given T and a small At, one must 
carefully choose the system size to be large enough such 
that the point (At, T) lies well above the crossover line. 

The model is exactly solvable in both the limits T — > 
and T — 5" oo. In the T — > limit, the thermal current J = 
Ei — E r is independent of L and the energy profile e(x) = 
\{Ei + E r ) is flat. In this limit Fourier law is violated as 
J is not proportional to local slope of e(x), which is now 
zero since e{x) is flat. The finite current J = E\ — E r is a 
consequence of the discontinuity of the energy profile at 
the boundaries. This discontinuity, and therefore a finite 
current, can never be obtained numerically for any finite 
At, however small. Numerical simulations, in fact, show 
that the current vanishes as J ~ At 2 for small At. 
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FIG. 11. (Color online) The crossover line T*(At) in the T- 
At plane for different L. In a finite system, this line separates 
the ballistic regime (shaded region below the curve), from the 
diffusive one. The data points are obtained numerically (see 
text) for L = 100, 200 (the broken line is only a guide to the 
eye). Following these trends, we schematically draw the same 
for very large L (solid line). Evidently, the apparent ballistic 
behaviour disappears for L — > oo and one obtains Fourier law 
for all T > 0, irrespective of the value of At. 



In the other solvable limit, i.e. when T — > oo, however, 
the energy profile e(x) = Ei + g £^' x is linear and one 
obtains a finite thermal conductivity k = 1. 

In conclusion, a thermodynamically large classical 
Heisenberg spin chain in one dimension obeys Fourier law 
for any non-zero temperature. However while studying 
thermal transport numerically for a finite system (though 
large) and a finite integration time step (though small) 
one must be careful in keeping the boundary temper- 
atures larger than the characteristic scale T*(At,L) to 
obtain the correct thermodynamic behaviour. Other- 
wise, for T<T* the system will show an apparent bal- 
listic behaviour which will eventually disappear in the 
L — > oo limit. This temperature dependent crossover 
from diffusive to ballistic behavior at small T is expected 
since T = is a critical point with a diverging correla- 
tion length. It will be quite interesting to study thermal 
transport in a system where one can set the boundary 
temperatures such that a singular point falls in the bulk. 
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